Effective single-particle order- N scheme for the dynamics of open non-interacting 

many-body systems 
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Quantum master equations are common tools to describe the dynamics of many-body systems 
open to an environment. Due to the interaction with the latter, even for the case of non- interacting 
electrons, the computational cost to solve these equations increases exponentially with the particle 
number. We propose a simple scheme, that allows to study the dynamics of N non-interacting 
electrons taking into account both dissipation effects and Fermi statistics, with a computational 
cost that scales linearly with N. Our method is based on a mapping of the many-body system 
to a specific set of effective single-particle systems. We provide detailed numerical results showing 
excellent agreement between the effective single-particle scheme and the exact many-body one, as 
obtained from studying the dynamics of two different systems. In the first, we study optically- 
induced currents in quantum rings at zero temperature, and in the second we study a linear chain 
coupled at its ends to two thermal baths with different (finite) temperatures. In addition, we give 
an analytical justification for our method, based on an exact averaging over the many-body states 
of the original master equations. 

PACS numbers: 03.65.Yz, 72.10.Bg 



I. INTRODUCTION 

Quantum systems that exchange energy with an envi- 
ronment have attracted a great deal of attention for many 
years 1,2 . The interest in these dissipative (open) quan- 
tum systems ranges from quantum computing and quan- 
tum information theory to biological physics 3 . Recent 
developments in the transport properties of nanoscale 
systems 4 raise new interest in these topics. For instance, 
the dissipative effects of the surrounding environment 
are key to understand the non-equilibrium properties 
of nanostructures and their approach to steady state 5 . 
However, the study of dissipative many-body quantum 
systems represents a major computational challenge. 

There are essentially two ways to approach this prob- 
lem. One consists in deriving equations of motion (mas- 
ter equations) for the reduced density matrix (DM) of the 
system of interest by integrating out the degrees of free- 
dom of the bath. 6 The further assumption of Markovian 
dynamics leads to different kinds of master equations for 
the DM 7 , perhaps the most popular being the Lindblad 
equation 8 which is often used in quantum optics 9,10 . The 
second approach is to use stochastic Schrodinger equa- 
tions 7,10 which are the stochastic unraveling of the mas- 
ter equations. If the Hamiltonian of the system does 
not depend on microscopic degrees of freedom, like the 
density or current density, both approaches describe the 
same physical properties. 11 

Irrespective of the chosen method, the solution of these 
equations is a formidable task which scales exponen- 
tially with the number of electrons. This is true even 
for a system of non-interacting electrons since the cor- 
relations induced by the bath make it impossible to ex- 
actly reduce the iV-particle equation of motion into N 
distinct single-particle equations of motion. It is the goal 
of this paper to discuss an ansatz which greatly simpli- 
fies this task for the dynamics of N non-interacting elec- 



trons in interaction with a bath. We focus on the DM 
approach but the conclusions are exactly the same for 
the stochastic Schrodinger equations. The latter, in fact, 
have found application in the recently developed stochas- 
tic time-dependent current-density functional theory (S- 
TDCDFT) 11 , an extension of time-dependent current- 
density functional theory to systems in dynamical inter- 
action with an environment. In S-TDCDFT the many- 
body interacting problem in the presence of the envi- 
ronment is mapped into an effective single-particle non- 
interacting problem in the presence of the same environ- 
ment. The ansatz we discuss in this work is thus of great 
use in the numerical solution of the equations of motion of 
S-TDCDFT, 11 and may therefore find application in dis- 
parate problems beyond the examples presented in this 
paper, where interactions are important. 



The outline of the paper is as follows. In Sec. II 
we describe in detail our proposed scheme. We define 
the master equation framework and our ansatz, along 
with the detailed structure of the resulting equations. 
In Sees. Ill and IV we give numerical examples of our 
scheme. We calculate currents induced by optical ex- 
citation in quantum ring structures in the presence of 
dissipation at T — (Sec. Ill) and consider steady state 
properties of a quantum system at finite temperatures 
(Sec. IV). We study systems which are small enough so 
that we can compare the results from our scheme with the 
full many-body calculation. We find excellent agreement 
between the two methods for a large range of parame- 
ters. In Sec. V we derive an analytical justification for 
our scheme. Starting from the exact many-body mas- 
ter equations we average over the many-body degrees of 
freedom and study the resulting (non-linear) equations. 
Sec. VI is devoted to a summary and outlook. 
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II. CALCULATION SCHEME 

Our goal is to study the dynamics of N electrons de- 
scribed by a non- interacting Hamiltonian TL — ■ Hj , 
while taking into account dissipation processes. To be 
more specific let us employ the following Lindblad-type 
master equation for the many-electron DM pu (h — 
1, e = I) 8 

Pm = -i[H, p M ] + £-Pm , (1) 

where [•] denotes the commutator, and C is the Lind- 
bladian superoperator, defined via a set V nn i of so-called 
"Lindblad operators" by 

CpM = {-\{ V L' V nn- , Pm} + V nn ,p M V^ , (2) 

with {•} the anti-commutator. The sums over n and n' 
(n 7^ n') are performed over all many-particle levels of 
the system and the F-operators are conveniently selected 
in the form V nn ' — ^/jnn'\^n}(^n'\, describing a transi- 
tion from the many-body state |* n ') m t° the state 
with the transition rate 7„„'. Although 7„„< are intro- 
duced phenomenologically here, these coefficients can be 
in principle derived from a microscopic theory. 

A common form for j nn / is described as follows 12 . At 
T = 0, dissipation drives the system towards its ground 
state, which we denote by the index n = 1. Therefore, it 
is reasonable to select 7„„< = for n > 1. Moreover, by 
assuming that the transition rate into the ground state is 
independent of n' , we may write 7i >n / = 7. This choice 
for the relaxation rates is a T = manifestation of de- 
tailed balance 7 , which we assume to hold for a Markovian 
ohmic bath in the long-time limit. In fact, there are other 
ways to choose the relaxation operators and still ensure 
detailed balance, and we have checked different options in 
our numerical calculations (Sec. Ill) and found no quali- 
tative change in our results. Therefore, we shall keep the 
above normalization hereon. 

For a system with M single-electron energy levels and 
N electrons, the solution of Eq. (1) generally requires the 
solution of (Cff + 2) x (Cff - l)/2 coupled differential 
equations, where Cff = M\/N\(M - N)\ and we have 
taken into account constrains of hermiticity and the unit 
trace of the density matrix. For the general case (exclud- 
ing, e.g., N = 1 or N = M), the problem thus scales 
exponentially with the number of particles 13 . 

Consider now an operator A = A? > a sum over 
single-particle operators. (This is not the most general 
form of operator but it encompasses most of the observ- 
ables of physical interest, like, e.g., the density or current 
density.) We make the following conjecture: the expec- 
tation value of A over a many-particle non-interacting 
electron state with dissipation can be approximated as a 
sum of single-electron expectation values of Aj over an 
ensemble of N single-electron systems with specifically 



selected single-electron dissipation operators, i.e. 

N 

TrApM^^TrA^K (3) 
i=i 

Here, p^ is a single-electron DM (effectively describing 
the j-th electron), each obeying its own Lindblad master 
equation 

= -i[H j ,pW]+ £«pW . (4) 

The choice of superoperators CS^ is dictated by two re- 
quirements: (i) for a time-independent Hamiltonian the 
dissipation processes should result in the Fermi-Dirac dis- 
tribution at long times, and (ii) the relaxation rate of 
many-electron states is 7. 

As we will demonstrate (numerically in Sec. Ill as well 
as analytically in Sec. V), these two requirements are 
met if one chooses a simple form for the T^-operators, 
which reflects the physical process at which the different 
electrons decay to consecutive single-particle levels (i.e., 
the i—th electron will decay to the i-th single-particle 
level, see Eq. (5)). Once a form for is chosen, one 
only needs to solve <~ N x M 2 equations, a reduction 
which enormously speeds up numerical calculations. 

The simplest choice for the Lindbladian superoperator 
which satisfies the above criteria is similar to the one in 
Eq. (2), with single-electron V operators of the following 
form: for the j-th electron we select at T = 

yi i = { Vl\3){k'\ , k'^k = j; k<k F 
kk ' y , otherwise ^ ' 

where \k) are now the single-particle states and ftp is the 
index of the Fermi level. In some sense, such a replace- 
ment of the many-body equation of motion by a set of 
auxiliary single-electron equations is similar to the intro- 
duction of a fictitious system of non-interacting electrons 
in density- functional theory. 11 

To summarize our scheme, it is constructed from the 
following steps: (i) given a non-interacting Hamiltonian, 
one constructs a set of Lindblad operators (following 
Eq. (2) and Eq. (5)), (ii) a set of single-particle density 
matrices p^ is defined, and corresponding master equa- 
tions [Eq. (4)] are solved, and, finally, (iii) any observ- 
able quantity (made of quadratic operators in the second 
quantization formalism) can be calculated using Eq. (3). 

III. NUMERICAL DEMONSTRATION: DRIVEN 
SYSTEM AT T — 

In order to test the conjecture (3), we have per- 
formed extensive numerical calculations considering a 
driven quantum system in a wide range of parameters. 
We found that for a system with non-degenerate levels 
Eq. (3) is almost perfectly satisfied. We believe that in 
systems with degenerate energy levels a deviation from 
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Eq. (3) is due to the intrinsic ambiguity of degenerate 
states. 

We consider a system of N tight-binding electrons on 
both a ring and a double ring of M sites in the presence of 
circularly-polarized electromagnetic radiation (see insets 
in Fig. 1). In order to lift the degeneracy, we place the 
system in a weak magnetic flux. The Hamiltonian of the 
system is given by 

H = -tJ2 (f^'^cla+x + fc.c.) +]T Uitycl* . (6) 

i i 

Here t is the hopping integral (we set \t\ = 1 through- 
out the calculation) and Ui{t) — —eE(t) • is a change 
of the potential energy of the i-th site (r^ is its po- 
sition) due to the external radiation. The magnetic 
field is taken into account via the usual Peierls substi- 
tution, with <p the magnetic flux through the ring, and 
<^o = h/e the flux quantum. The electric field is written 
as E(t) = Eq cos(iot)-k±E sin(wi)y, where E and u> are 
the electric field amplitude and frequency, x and y are 
unit vectors in the x and y directions (in the ring plane) , 
and ± corresponds to a a± circular polarization. 

It is known that in the ring topology a circularly- 
polarized radiation creates a current in the ring 14 ' 15 . We 
calculate the expectation value of the current operator 
through a specific bond, J = ¥(c]ct+i — h.c), using both 
the exact many-body DM, and a set of single-electron 
density matrices calculated as described above. 16 In both 
schemes we start by diagonalizing the tight-binding part 
of the Hamiltonian. In the many-body (exact) scheme, 
we then write the time-dependent potential and the Lind- 
blad operators in their full many-body form and solve the 
time-dependent set of equations for the many-body DM. 
For the single-particle scheme, we solve a set of N single- 
particle Lindblad equations (of size M x M), each with 
its own set of relaxation operators D^'. The current 
is then calculated as a function of time using the LHS 
(many-body form) and the RHS (single-particle form) of 
Eq. (3). The calculations were made for a wide range 
of system parameters, displaying excellent agreement be- 
tween the two schemes. 



A. Ground-state initial conditions 

Fig. 1 shows the current calculated by the two meth- 
ods through a bond connecting two adjacent sites of a 
10-site ring containing 3 electrons. We see that the cur- 
rent through the bond oscillates in agreement with a pre- 
vious study 15 . Most importantly, in the context of the 
present investigation, the current values hardly differ be- 
tween the two schemes. The average deviation of the 
two currents is less then 1.5% (the maximum deviation 
is « 6.5%). This difference rapidly disappears with de- 
creasing Eq. This is seen from comparing Fig. 2(a) and 
2(b) where the current excited in a double ring struc- 
ture is plotted for two different values of the electric field 
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FIG. 1: (Color online) Current between two sites of a ring as 
a function of time calculated by the exact many-body and ap- 
proximate single-electron approaches. Inset shows the system 
geometry. This calculation has been done with the following 
set of parameters: N = 3, M — 10, eEoa = 0.1, u) = 0.8, 
a = 1, 7 = 0.1, a = 0.1415nm and B = 10T. a is a bond 
length. The magnetic field corresponds to a flux through the 
ring of (j)/(j> « 1.66 x 10" 4 . 

amplitude, eEoa = 0.1 and eEoa — 0.01, respectively. 
One clearly sees that the discrepancies between the two 
methods (marked in a gray circle in 2(a)) diminish as the 
excitation field decreases. The results presented in Figs. 
1 and 2 were obtained assuming that at time t = the 
system is in its ground state. 

B. Non-equilibrium initial conditions 

Next, we have tested the applicability of our approach 
to highly excited states. In Fig. 3 we plot the current 
gnerated in a 10-site ring containing 3 electrons. The dif- 
ference with previously discussed calculations is that now 
we assume that in the initial moment of time the system 
is in its highest energy state. For a single electron, there 
are 10 energy states in the 10-site ring. We made calcu- 
lations considering different relaxation schemes. Indeed, 
there is an arbitrariness in the relaxation state assign- 
ment (e.g., V operators for the electron which is initially 
in the 10-th state - highest energy state - can be se- 
lected to describe its relaxation into the first, second or 
third lowest energy state). Fig. 3 displays a very good 
agreement of the many-body calculation compared to the 
results obtained using our single-electron approach with 
two different relaxation schemes. Importantly, since the 
rates at which electrons relax into their ground states are 
the same, the two relaxation schemes lead to the same 
current, showing the insensitivity of our general scheme 
to the details in the initial state de-population. Also, in 
the long-time limit the current is independent of the ini- 
tial conditions chosen (cf. the current in Fig. 3 with the 
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FIG. 2: (Color online) Current excited in a double ring cal- 
culated by the exact many-body and approximate single- 
electron approaches. The electric field amplitude is (a) 
eEoa = 0.1 and (b) eEoa = 0.01. All other parameters are the 
same as in Fig. 1. The system geometry is shown in the inset 
of (a). The discrepancies between the two methods (such as 
marked in a gray circle in (a)) diminish as the excitation field 
decreases. 



current in Fig. 1 at t > 40). 



C. Precision of the simplified scheme 

In order to study the precision of the single-electron 
scheme, we calculate the current excited in a 6-site ring 
containing 3 electrons. Our main observation is that the 
simplified scheme provides a very good precision for the 
whole range of parameters used in the calculations. We 
have found that a slightly better precision is obtained at 
weak and strong electric fields. This particular observa- 
tion is clearly seen in Fig. 4 where we plot the ratio of 
the RMS of current differences calculated as 



Ajrms f \J (jmb Jse) 



dt 



(7) 
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FIG. 3: (Color online) Comparison of many-body calcula- 
tions with those obtained with the simplified approach in the 
case of non-equilibrium initial conditions. For single electron 
calculations, we used two different relaxation schemes shown 
as insets. This plot was obtained for a 10-site quantum ring 
using the same parameter values as in Fig. 1. 
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FIG. 4: (Color online) RMS of the difference of currents cal- 
culated by many-body and single-electron schemes divided by 
the maximum current amplitude within the calculation time 
as a function of the electric field amplitude. This plot was 
obtained for N = 3, M = 6, a = 1, a = 0.1415nm, B = 10T 
and r = 100. The other calculation parameters are shown in 
the figure. 



to Aj max = jmb x ~ Jmb"- Here > T 58 a sampling pe- 
riod, jmbtse) IS the current calculated using many-body 



max(min) 



is the maximum 



(single-electron) scheme and j 
(minimum) value of current calculated within the time 
interval [0, r]. A better agreement at weak fields can be 
related to the fact that in this situation only the low- 
energy states become occupied and the relaxation oper- 
ators in the many-electron and single-electron schemes 
are the same (see Sec. V for more arguments). At strong 
fields, the better agreement is due to the fact that the 
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electric field term is dominant in the equations of mo- 
tion. Fig. 4 also demonstrates that the single-electron 
scheme precision slightly depends on simulation param- 
eters and is a better approximation when dissipation is 
weaker. 

Fig. 5 presents selected results of our calcula- 
tions showing agreement between many-body and single- 
electron calculations at several values of the electric field 
amplitude. The interesting feature of these results is that 
at weak driving fields the single-electron scheme precision 
is better at longer times (t > 60 in Fig. 5(a)), at inter- 
mediate fields the scheme precision is better in the initial 
time interval (t < 20 in Fig. 5(b)) and at strong fields 
the precision is better again at longer times (Fig. 5(c)). 



IV. NUMERICAL DEMONSTRATION: STEADY 
STATE AT FINITE TEMPERATURES 

In the second numerical example, we study a non- 
equilibrium system at finite temperatures. The system 
of interest is a linear metallic chain, connected at its two 
ends to two thermal baths at different temperatures, Tl 
and Tr, corresponding to the left and right temperatures 
(see inset of Fig. 6(b)). 

The Hamiltonian of the system is given by TL = 

~tY!,(i,j)eL,R,d (4 c i +h.c)j (t is the hopping integral, 
which serves as the energy scale, and we have chosen 
t = 1). The master equation now takes the form 



-i[H, p] + C L [ P ] + C R [p] 



(8) 



where Cl(r) describes relaxation processes due to the 
contact between the left (right) lead with its respective 
bath at temperature ?l(_r)- The ^-operators are given 
by 



(L,R) 



kk' 



^ R) f { n L ' R \e k )\k)(k'\ 



(9) 



where f£' R) (e k ) = 1/ (exp ( k £ ^ R ) + l) are the Fermi 
distributions of the left and right leads, with p the chem- 
ical potential. The coefficients 



(L,R) 
Ikk' 



\i>k{r)~fO^*k>( r )\r=r L {rR) 



(10) 



describe the overlap between the single-particle states 
\k) and \k') over the point of contact rum between the 
left (right) baths and the corresponding junction leads. 
The constant 70 describes the strength of interactions 
between the bath and electrons. The form (10) can be 
derived from first principles by tracing out the bath de- 
grees of freedom, with the latter formed by a dense spec- 
trum of boson excitations (e.g., phonons), which inter- 
act locally with electrons at the edges of the system. 
Physically, it corresponds to the experimental situation 
in which the left (right) bath induces energy relaxation 
only between states which reside predominantly on the 
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FIG. 5: (Color online) Current excited in a 6-site ring cal- 
culated by different approaches as indicated. The calculation 
parameters are the same as in Fig. 4. The electric field ampli- 
tude is eEoa — 0.1 (a), 1.7 (b) and 8 (c). The insets show the 
absolute value of many-body (Jmt>) and single-electron (j S e) 
currents difference as a function of time. 



left (right) edge of the junction, where the bath is in 
contact. The operators (9) guarantee that the system 
evolves to a global equilibrium if Tl = Tr. For Tl =^ Tr 
this system is inherently out of equilibrium, and reaches a 
steady state which may have, for instance, a non-uniform 
electron density 17 , and is thus relevant for experiments 
of thermo-power measurements in nano-systems 18 . We 
point out that the above model also relaxes the con- 
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FIG. 6: (Color online) Occupation of the different single- 
particle energy levels as a function of time for the two cal- 
culation schemes, the full many-body calculation (solid lines) 
and the approximate scheme (dashed lines). Initial condition 
are either (a) the ground state or (b) a uniformly-occupied 
state. The chain length is L — 12, with the parameters g = 1, 
70 = 0.01, T L = 0.1 and T R = 0.4. 
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FIG. 7: (Color online) Difference in the local density An (av- 
eraged over the entire chain) as a function of system length 
at steady state. The numerical parameters are the same 
as in Fig. 6. Inset : local density along a L = 16 chain, 
calculated using the exact scheme (points) and approximate 
scheme (solid line). 



for the same parameters as in Fig. 6 for different chain 
lengths L = 5,6, 16. We find that as the system be- 
comes larger the approximation improves (the relative 
deviation for the larger systems is less than 3%). The 
reason for the improvement of the approximation with 
increasing length stems from the fact that as the sys- 
tem becomes larger, the single-particle occupations of 
the many-body system become closer and closer to a true 
broadened Fermi distribution. In the inset of Fig. 7 we 
plot the local density along a L — 16 chain, calculated 
using the exact scheme (points) and approximate scheme 
(solid line), showing the excellent agreement between the 
two. 



straint of Sec. Ill that there is a single relaxation rate 
for all relaxation processes. 

In Fig. 6 we plot the occupation of the different single- 
particle energy levels as a function of time for the two 
calculation schemes, the full many-body (solid lines) and 
the approximate scheme (dashed lines). The chain length 
is L = 12, with the parameters g = 1, 70 = 0.01, 
Tl = 0.1 and Tr = 0.4, and it is occupied by two elec- 
trons. We have plotted the dynamics starting from either 
the ground state (Fig. 6(a)) or a uniform state, where all 
energy levels are equally occupied (Fig. 6(a)). As seen, 
starting from the ground state (Fig. 6(a)) there is ex- 
cellent agreement between the two schemes both in the 
transient dynamics and in the steady state. On the other 
hand, if we start from an excited state (Fig. 6(b)) then 
the transient dynamics exhibit slight differences between 
the exact and approximate scheme. The steady state is, 
naturally, the same with either initial conditions. Sim- 
ilar calculations with different parameters have yielded 
similar results. 

In order to study the accuracy of the approxima- 
tion also in the present example, we calculate the dif- 
ference in the local density between the two schemes, 
Ani = \nt rnb - rii >S p\, at steady state. Here, n iiTO i,( sp ) is 
the local density (n, = J2 k \^pk{i)\ 2 pkk) at the i-th site, 
calculated with the many-body (single-particle) scheme. 
In Fig. 7 we plot An (averaged over the entire chain), 



V. ANALYTIC JUSTIFICATION 

We now provide an analytical argument for the validity 
of our ansatz , which is summarized in Eqs. (3-5) for 
T = 0. In order to do so we start from the definition of an 
auxiliary single-particle density matrix (SPDM) from the 
many-body one. We then evaluate its equation of motion 
by summing up the many-body degrees of freedom, and 
study the structure of the equations. We in fact find 
that this SPDM can be approximately written as sum 
of single-particle density matrices obeying equations of 
motion with specific bath operators, thus validating our 
ansatz. We do this for finite temperatures, and show 
that the result leads to the T = form for the relaxation 
operators used in the numerical calculations. 

Let us define the following SPDM 

p{t)=Y,Pkk-{t)\k){k'\ . (11) 

kk' 

The matrix elements are derived from the many-body 
DM by 

Pkk' = Tr (clck'PM) ■ (12) 
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We show below that p{t) can be approximated as 

p(t)^J2p U) (t), (13) 

3 

where p^ are the single-particle density matrices enter- 
ing Eq. (3). 

The time evolution of the SPDM is determined by 

Pkk' = ^;Tr (c\c kl pM^ = Tr (c\c k '{-i\H, p M ] + £/°m)) ■ 

, (14) 

One can now perform the trace exactly using Wick's 
theorem. The relaxation operators V nn > defined below 
Eq. (2) generally involve up to M creation and M an- 
nihilation operators. Therefore, it is not practical to 
use them in analytical calculations. We instead con- 
sider V operators of a commonly used 19 simplified form 
l , 

Vkk' = (7fcfc') 2 c\ck' ,k^k. It is clear that when exci- 
tation of the system is weak, and highly excited states are 
almost unpopulated, the physical effect caused by both 
operators is nearly the same. Note, however, that tak- 
ing this form for the ^-operators (which excludes direct 
relaxation of highly-excited states into the ground state) 
does not lead to a reduction in the number of equations 
needed to be solved, since the equations remain fully cou- 
pled (put it differently, the Lindbladian operator cannot 
be subdivided into blocks). 

A. Diagonal elements 

We start by deriving equations of motion for the di- 
agonal elements of the SPDM. For the sake of simplic- 
ity, let us assume that the system Hamiltonian is time- 
independent and diagonalizcd. Then, it is easy to find 
that the equations of motion for the diagonal elements of 
the SPDM are 20 

Pkk — — ^ X/ Tfcfc'Pfefe + ^ X] Ik'kPk'k 1 + 
k'^k k'^k 

+ ^Pkk ^2 (jk'k ~ lkk' )Pk'k' ■ (15) 
k'jtk 

Let us examine Eq. (15) by making two assump- 
tions: (i) the coefficients are only a function of the first 
index, i.e., ^ k k> — 7fc', and (ii) the third (non-linear) part 
on the RHS of Eq. (15) is negligible and is set to zero. 
Within these assumptions, and noting that by definition 
Sfeli Pkk = N, one obtains the equation 

pkk = -Zp k k+Jk(N - pkk), (16) 

where Z = J2 k ,^ k 7k'- Solving this equation yields 

Pkkit) = (pkk(0) - e-<*^>« + (17) 

V Ik + z J lk + z 

For a Fermi system, the long-time limit of the SPDM 
should be p kk (t — » oo) = /c(efc), where /c(efe) = 1/(1 + 



cxp((efe — p)/ksT) is the Fermi-Dirac distribution. It 
follows directly that in order to satisfy this long-time 
limit, the coefficients must be chosen such that j k = 
7/u(efe)- 

We now turn back to the third, non-linear part in the 
RHS of Eq. (15). Keeping in mind the definition for j k , 
this part now reads jp kk Efc'/fc(.M £ fc) - /£>(efc'))Pfc'fc'- 
In the long-time limit, as pkk approach their equilibrium 
values, and at zero temperature, one can consider two 
possibilities. In the first, both k and k' lie below or above 
the Fermi surface. In this case, /z)(efc) — fo{^k') ~ and 
the non-linear part vanishes. If, on the other hand, cither 
k or k' lie below the Fermi surface and the other above it, 
then indeed /n(efe) — fo(^k') 7^ 0. However, in that case 
either p kk w or p k iy rts 0. Thus, in the low temperature 
long-time limit, the third term on the RHS of Eq. (15) is 
negligible, which means that our assumption (ii) above 
is justified. 

Extending this conclusion to finite temperatures and 
to all times, we end up with a simple equation for the 
diagonal elements of the SPDM, 

Pkk = ~7 ^ fD{ek')Pkk +1^2 fD( e k)pk'k'- (18) 

k'yik k'^k 

Simple algebra reveals that these equations are equal to 
those obtained from applying the Lindbladian operator, 
Eq. (2), to the SPDM, with the ^-operators having the 
form 

V kk ' = VjfD(ek)\k)(k'\ , (19) 

which is a particular case of the operators (9), thus jus- 
tifying their structure. We thus propose that the SPDM 
evolves according to Eq. (1) and (2), with the Lindblad 
operator given in terms of Eq. (19). 

The equations of diagonal SPDM elements can be de- 
rived differently. Since p(t) ~ ^2jP^(t), using Eq. (4) 
with the single-electron ^-operators in the form 

Kk> = Ml - Skk>)Vlh^k)\3)(k'l (20) 

we can obtain a set of equations which is the same as 
Eq. (15). This demonstration clearly shows a similarity 
of our single-electron and many-body approaches. Note, 
that the definition (20) coincides with Eq. (5) at T = 0. 
Moreover, while there is no a priori justification for ne- 
glecting the non-linear terms, the numerical calculations 
of the previous sections show that it is an excellent ap- 
proximation for non-interacting systems. 

Let us also point out that the equations for the diag- 
onal and off-diagonal parts of the SPDM are completely 
decoupled (this result is exact). Therefore, if one is in- 
terested in the time-dependent expectation value of an 
operator that commutes with the Hamiltonian, or only 
in the steady-state (where the off-diagonal elements van- 
ish) our ansatz reduces the computational effort to a sin- 
gle M x M equation for the diagonal elements of the 
SPDM. 
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B. Off-diagonal elements 

The off-diagonal elements of the density matrix are 
needed to calculate, e.g., local currents or densities in 
a non-equilibrium situation of an excited system (as in 
the numerical examples of Sec. III). As stated above, 
if only the diagonal elements are of interest, SPDM cal- 
culations with the ^-operators in their especially simple 
form (Eq. 19) can be used. If the off-diagonal elements 
are important, calculations using single-electron matrices 
p 3 with relaxation operators given by Eq. 20 have to be 
performed. 

In order to understand why single-electron calculations 
are needed (or why SPDM does not provide the best 
results in all cases), we study the equation of motion for 
the off-diagonal elements of the exact Lindblad operator. 
Using Eq. (14) and Vkk 1 operators defined below Eq. (14) 
one finds 

(jCp)fefe' = -- ^2 {lk"k> + 7fc"fc)(l - Pk"k")pkk' - 

— - ^2 ilkk" + lk'k")Pkk' Pk"k" ■ (21) 
k"^k,k' 

Again we make the substitution "fkk 1 — lfD{^k')} an< i 
consider for simplicity the system at zero temperature. 
By assuming pi~k ~ /i)(e/c) we find that the first element 
of the LHS in Eq. (21) is negligible, and one is left with 



kk' 



1 y 

2 ^ 

k"^k,k' 



Pk"k"Pkk' = —-^(N—pkk — pk'k')pkk' 

(22) 

If one uses SPDM calculations to study the off-diagonal 
elements, then one finds that (Cp)kk' does not depend 
on pkk, Pk'k' at all. However, within the single-electron 
scheme this separation can not be made, and the dynam- 
ics of the off-diagonal elements are better captured. This 
can be seen in the numerical example by comparing the 
exact many-body calculation with the approximate cal- 
culation using both Eq. (13) and the SPDM Eq. (11). 
This is shown in Fig. 8, where a comparison between the 
three methods is shown. As seen in the figure, the agree- 
ment between all schemes is good in general, with sub- 
stantial differences arising only at the maxima and min- 
ima of the current. At these points, the single-particle 
scheme [Eq. (13)] is closer to the many-body calculation 
than the SPDM method [Eq. (11)]. 



VI. SUMMARY 



(Eq. (4)) where both Fermi statistics and dissipation are 
1 1 ' 1 1 1 1 1 
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FIG. 8: (Color online) Current excited in a 10-site ring cal- 
culated by three different methods as indicated in the figure. 
This plot was obtained using the same parameter values as in 
Fig. 1 except u = 1. 



taken into account via a specific form of relaxation oper- 
ators (Eq. (5) at T = 0; Eq. (20) for T ^ 0). We have 
numerically demonstrated that the proposed method is 
in excellent agreement with the exact many-body calcu- 
lation by studying two example. The first example is a 
system of tight-binding rings at zero temperature, driven 
out of equilibrium by external radiation. The second ex- 
ample is a linear chain connected at its end to two heat 
baths held at different temperatures. 

Since, even for non-interacting electrons the inclusion 
of the Pauli exclusion principle is nontrivial for open 
quantum systems 21 , we believe our scheme can be used 
in systems where interactions play a relatively minor 
role such as in graphene 22 , quantum point contacts 23 , 
etc. Nevertheless, while the above examples did not 
include electron-electron interactions, the latter may 
be included within the framework of stochastic time- 
dependent current-density functional theory 11 , where the 
interacting many-body problem in the presence of envi- 
ronments is mapped into an effective single-particle prob- 
lem in the presence of the same environments. Our ansatz 
thus provides a good starting point to solve the corre- 
sponding equations of motion with a computational cost 
that scales only linearly with the number of particles. 
Such a project is currently underway. 



We have proposed an order- N scheme to investigate 
the dynamics of N non-interacting electrons coupled to 
one or more baths, and justified it analytically by examin- 
ing and tracing the full many-body calculation. The main 
idea is to reduce the equation of motion for the many- 
body system to a set of effective single-electron equations 
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